1
2
- 2.1 OPTICS
- 2.2 DBSCAN
3
4
- 4.1 Ventajas k-means y DBSCAN
- 4.2 Medidas para mitigar las desventajas

El estudio se realizará en base al juego de datos Hawks presente en el paquete R Stat2Data.

Los estudiantes y el profesorado del Cornell College en Mount Vernon, Iowa, recogieron datos durante muchos años en el mirador de halcones del lago MacBride, cerca de Iowa City, en el estado de Iowa. El conjunto de datos que analizamos aquí es un subconjunto del conjunto de datos original, utilizando sólo aquellas especies para las que había más de 10 observaciones. Los datos se recogieron en muestras aleatorias de tres especies diferentes de halcones: Colirrojo, Gavilán y Halcón de Cooper.

Hemos seleccionado este juego de datos por su parecido con el juego de datos penguins y por su potencial a la hora de aplicarle algoritmos de minería de datos no supervisados. Las variables numéricas en las que os basaréis son: Wing, Weight, Culmen, Hallux

data("Hawks")
summary(Hawks)

##      Month             Day             Year       CaptureTime   ReleaseTime 
##  Min.   : 8.000   Min.   : 1.00   Min.   :1992   11:35  : 14          :842  
##  1st Qu.: 9.000   1st Qu.: 9.00   1st Qu.:1995   13:30  : 14   11:00  :  2  
##  Median :10.000   Median :16.00   Median :1999   11:45  : 13   11:35  :  2  
##  Mean   : 9.843   Mean   :15.74   Mean   :1998   12:10  : 13   12:05  :  2  
##  3rd Qu.:10.000   3rd Qu.:23.00   3rd Qu.:2001   14:00  : 13   12:50  :  2  
##  Max.   :11.000   Max.   :31.00   Max.   :2003   13:05  : 12   13:32  :  2  
##                                                  (Other):829   (Other): 56  
##       BandNumber  Species  Age     Sex          Wing           Weight      
##            :  2   CH: 70   A:224    :576   Min.   : 37.2   Min.   :  56.0  
##  1142-09240:  1   RT:577   I:684   F:174   1st Qu.:202.0   1st Qu.: 185.0  
##  1142-09241:  1   SS:261           M:158   Median :370.0   Median : 970.0  
##  1142-09242:  1                            Mean   :315.6   Mean   : 772.1  
##  1142-18229:  1                            3rd Qu.:390.0   3rd Qu.:1120.0  
##  1142-19209:  1                            Max.   :480.0   Max.   :2030.0  
##  (Other)   :901                            NA's   :1       NA's   :10      
##      Culmen         Hallux            Tail        StandardTail  
##  Min.   : 8.6   Min.   :  9.50   Min.   :119.0   Min.   :115.0  
##  1st Qu.:12.8   1st Qu.: 15.10   1st Qu.:160.0   1st Qu.:162.0  
##  Median :25.5   Median : 29.40   Median :214.0   Median :215.0  
##  Mean   :21.8   Mean   : 26.41   Mean   :198.8   Mean   :199.2  
##  3rd Qu.:27.3   3rd Qu.: 31.40   3rd Qu.:225.0   3rd Qu.:226.0  
##  Max.   :39.2   Max.   :341.40   Max.   :288.0   Max.   :335.0  
##  NA's   :7      NA's   :6                        NA's   :337    
##      Tarsus        WingPitFat        KeelFat           Crop       
##  Min.   :24.70   Min.   :0.0000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:55.60   1st Qu.:0.0000   1st Qu.:2.000   1st Qu.:0.0000  
##  Median :79.30   Median :1.0000   Median :2.000   Median :0.0000  
##  Mean   :71.95   Mean   :0.7922   Mean   :2.184   Mean   :0.2345  
##  3rd Qu.:87.00   3rd Qu.:1.0000   3rd Qu.:3.000   3rd Qu.:0.2500  
##  Max.   :94.00   Max.   :3.0000   Max.   :4.000   Max.   :5.0000  
##  NA's   :833     NA's   :831      NA's   :341     NA's   :343

1

Presento el juego de datos, nombre y significado de cada columna, así como las distribuciones de sus valores.

Adicionalmente realizo un estudio tipo EDA (exploratory data analysis) similar al de los ejemplos 1.1 y 1.2 ( k-means )

1.1 Exploración de los datos

Miramos la estructura de los datos para ver el nombre de cada columna y la definición de sus datos. Las variables numéricas en las que nos basaremos son:

Wing, Medición del ala.
Weight, Medición del peso.
Culmen, Medida de culmen de cada halcón.
Hallux, Medida del primer dedo de la pata en cuanto a orden.
Species, Etiqueta de la especia de halcón.

str(Hawks)

## 'data.frame':    908 obs. of  19 variables:
##  $ Month       : int  9 9 9 9 9 9 9 9 9 9 ...
##  $ Day         : int  19 22 23 23 27 28 28 29 29 30 ...
##  $ Year        : int  1992 1992 1992 1992 1992 1992 1992 1992 1992 1992 ...
##  $ CaptureTime : Factor w/ 308 levels " ","1:15","1:31",..: 181 25 138 42 62 71 181 88 261 192 ...
##  $ ReleaseTime : Factor w/ 60 levels ""," ","10:20",..: 1 2 2 2 2 2 2 2 2 2 ...
##  $ BandNumber  : Factor w/ 907 levels " ","1142-09240",..: 856 857 858 809 437 280 859 860 861 281 ...
##  $ Species     : Factor w/ 3 levels "CH","RT","SS": 2 2 2 1 3 2 2 2 2 2 ...
##  $ Age         : Factor w/ 2 levels "A","I": 2 2 2 2 2 2 2 1 1 2 ...
##  $ Sex         : Factor w/ 3 levels "","F","M": 1 1 1 2 2 1 1 1 1 1 ...
##  $ Wing        : num  385 376 381 265 205 412 370 375 412 405 ...
##  $ Weight      : int  920 930 990 470 170 1090 960 855 1210 1120 ...
##  $ Culmen      : num  25.7 NA 26.7 18.7 12.5 28.5 25.3 27.2 29.3 26 ...
##  $ Hallux      : num  30.1 NA 31.3 23.5 14.3 32.2 30.1 30 31.3 30.2 ...
##  $ Tail        : int  219 221 235 220 157 230 212 243 210 238 ...
##  $ StandardTail: int  NA NA NA NA NA NA NA NA NA NA ...
##  $ Tarsus      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ WingPitFat  : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ KeelFat     : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Crop        : num  NA NA NA NA NA NA NA NA NA NA ...

Nos encontramos con 908 observaciones con 19 variables. Vemos datos faltantes en alguna variable de interés.

1.2 Limpieza

#Revisar valores non-NA en cada columna
colSums(is.na(Hawks))

##        Month          Day         Year  CaptureTime  ReleaseTime   BandNumber 
##            0            0            0            0            0            0 
##      Species          Age          Sex         Wing       Weight       Culmen 
##            0            0            0            1           10            7 
##       Hallux         Tail StandardTail       Tarsus   WingPitFat      KeelFat 
##            6            0          337          833          831          341 
##         Crop 
##          343

Encontramos datos faltantes en todas las variables de interés:

#Revisar valores vacíos
void_data<-which(colSums(Hawks=="")!=0)
names(void_data)

## [1] "ReleaseTime" "Sex"

Existen valores en blanco en variables que no nos influye para nuestro análisis.

Seleccionamos las variables de interés

data <- Hawks[,c("Species","Wing", "Weight", "Culmen", "Hallux")]
data[1:5, 1:5]

##   Species Wing Weight Culmen Hallux
## 1      RT  385    920   25.7   30.1
## 2      RT  376    930     NA     NA
## 3      RT  381    990   26.7   31.3
## 4      CH  265    470   18.7   23.5
## 5      SS  205    170   12.5   14.3

1.3 Estudio por variable de interés

1.3.1 Wing

# Resumen estadístico de Wing
summary(data$Wing)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    37.2   202.0   370.0   315.6   390.0   480.0       1

Tenemos un valor nulo.

# Generamos histograma de Wing
hist(data$Wing)

La distribución del atributo Wing no es normal, es bimodal con mayor frecuencias de medidas de ala entre 350-400 y 150-200.

# Generamos boxplot por etiqueta
boxplot(data$Wing, main =colnames(data$Wing), 
        xlab = "Especies", ylab = paste("Mediciones",colnames(data$Wing)))

No encontramos outliers en el atributo Wing. Ahora paso a imputar la mediana a los valores nulos.

# Reemplazo valores nulos por la mediana
data$Wing[is.na(data$Wing)] <- median(data$Wing, na.rm = TRUE)
sum(is.na(data$Wing))

## [1] 0

Ahora tenemos el atributo Wing sin ningún valor nulo.

1.3.2 Weight

# Resumen estadístico de Wing
summary(data$Weight)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    56.0   185.0   970.0   772.1  1120.0  2030.0      10

Encontramos 10 valores nulos.

# Generamos histograma de Wing
hist(data$Weight)

La distribución del atributo Wing no es normal, es bimodal con mayor frecuencias de peso entre 0-200 y 1000-1200.

# Generamos boxplot por etiqueta
boxplot(data$Weight, main =colnames(data$Weight), 
        xlab = "Especies", ylab = paste("Mediciones",colnames(data$Weight)))

No encontramos outliers en el atributo Weight. Ahora paso a imputar la mediana a los valores nulos.

# Reemplazo valores nulos por la mediana
data$Weight[is.na(data$Weight)] <- median(data$Weight, na.rm = TRUE)
sum(is.na(data$Weight))

## [1] 0

Ahora tenemos el atributo Weight sin ningún valor nulo.

1.3.3 Culmen

# Resumen estadístico de Culmen
summary(data$Culmen)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     8.6    12.8    25.5    21.8    27.3    39.2       7

Encontramos 7 valores nulos.

# Generamos histograma de Culmen
hist(data$Culmen)

La distribución del atributo Culmen no es normal, es bimodal con mayor frecuencias de medida entre 24-27.

# Generamos boxplot por etiqueta
boxplot(data$Culmen, main =colnames(data$Culmen), 
        xlab = "Especies", ylab = paste("Mediciones",colnames(data$Culmen)))

No encontramos outliers en el atributo Culmen. Ahora paso a imputar la mediana a los valores nulos.

# Reemplazo valores nulos por la mediana
data$Culmen[is.na(data$Culmen)] <- median(data$Culmen, na.rm = TRUE)
sum(is.na(data$Culmen))

## [1] 0

Ahora tenemos el atributo Culmen sin ningún valor nulo.

1.3.4 Hallux

# Resumen estadístico de Hallux
summary(data$Hallux)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    9.50   15.10   29.40   26.41   31.40  341.40       6

Enocntramos 6 valores nulos.

# Generamos histograma de Hallux
hist(data$Hallux)

La distribución del atributo Hallux no es normal, tiene cola a derecha.

# Generamos boxplot por etiqueta
boxplot(data$Hallux, main =colnames(data$Hallux), 
        xlab = "Mediciones", ylab = paste("Hallux",colnames(data$Hallux)), horizontal = TRUE)

Apreciamos outliers en la variable Hallux. Eliminamos todos los valores superiores a 70 son todos los outliers encontrados en el atributo.

data <- data[-which(data$Hallux > 70),]

boxplot(data$Hallux, main =colnames(data$Hallux),
xlab = "Mediciones", ylab = paste("Hallux",colnames(data$Hallux)), horizontal = TRUE)

Hemos eliminado los outliers, paso a imputar la mediana a los valores nulos.

# Reemplazo valores nulos por la mediana
data$Hallux[is.na(data$Hallux)] <- median(data$Hallux, na.rm = TRUE)

summary(data)

##  Species       Wing           Weight           Culmen          Hallux     
##  CH: 70   Min.   : 37.2   Min.   :  56.0   Min.   : 8.60   Min.   : 9.50  
##  RT:574   1st Qu.:203.0   1st Qu.: 188.0   1st Qu.:12.80   1st Qu.:15.10  
##  SS:257   Median :370.0   Median : 970.0   Median :25.50   Median :29.40  
##           Mean   :316.0   Mean   : 775.5   Mean   :21.85   Mean   :25.27  
##           3rd Qu.:390.0   3rd Qu.:1120.0   3rd Qu.:27.30   3rd Qu.:31.40  
##           Max.   :480.0   Max.   :2030.0   Max.   :39.20   Max.   :54.50

Creo matriz de correlación, ofrece una forma rápida de comprender la fuerza de las relaciones lineales que existen entre las variables en el conjunto de datos.

Estudio de la correlación

#create correlation matrix
cor(data[,2:5])

##             Wing    Weight    Culmen    Hallux
## Wing   1.0000000 0.9245828 0.9506569 0.9165801
## Weight 0.9245828 1.0000000 0.9456756 0.9126963
## Culmen 0.9506569 0.9456756 1.0000000 0.9455982
## Hallux 0.9165801 0.9126963 0.9455982 1.0000000

Las correlaciones son muy altas, dominarán el cálculo de distancia y el resultado de agrupación será más dependiente de ellos, lo que no se desea. Si nuestras variables fueran totalmente independientes no habría ningún problema. Sin embargo, si tienen algún tipo de correlación, una influye sobre la otra y esta influencia no queda bien reflejada si usamos la distancia estadística. Para corregir la distorsión provocada por la correlacionados, la agrupación basada en la distancia de Mahalanobis sería la más apropiada.

Estudio de la covarianza

cov(data[,2:5])

##              Wing     Weight     Culmen     Hallux
## Wing    9030.6679  40333.085  654.42216  721.46504
## Weight 40333.0852 210722.719 3144.64854 3470.29892
## Culmen   654.4222   3144.649   52.47452   56.73693
## Hallux   721.4650   3470.299   56.73693   68.60722

Los valores a lo largo de las diagonales de la matriz son simplemente las varianzas de cada variable. Los otros valores de la matriz representan las covarianzas entre los distintas variables. Las cuatro variables tienen una covarianza positiva. Un número positivo para la covarianza indica que dos variables tienden a aumentar o disminuir en tándem. Por ejemplo, Weight y Wing tienen una covarianza positiva (40333.085), lo que indica que los halcones que obtienen una medición alta en peso también tienden a obtener una medición alta en tamaño de las alas. Por el contrario, los hacones que con poco peso también tienden a obtener medidas bajas de las alas.

Estandarización de los datos

Ahora escalamos datos para que las variables tengan una media de 0 y una desviación estándar de 1. De esta evitamos diferencias de escala, que alguna variable tenga mayor peso que otra, pero sin influir en las diferencias entre variables o sea que la correlación seguirá siendo alta.

# normalizamos columnas objetivo
data.norm <- scale(data[2:5])
head(data.norm)

##         Wing     Weight     Culmen     Hallux
## 1  0.7262840  0.3147486  0.5308972  0.5826527
## 2  0.6315769  0.3365329  0.5032879  0.4981418
## 3  0.6841920  0.4672389  0.6689439  0.7275287
## 4 -0.5364774 -0.6655465 -0.4354291 -0.2141650
## 5 -1.1678581 -1.3190766 -1.2913181 -1.3248807
## 6  1.0104054  0.6850823  0.9174278  0.8361857

# comprovamos media y var
summary(data.norm)

##       Wing             Weight            Culmen            Hallux       
##  Min.   :-2.9336   Min.   :-1.5674   Min.   :-1.8297   Min.   :-1.9044  
##  1st Qu.:-1.1889   1st Qu.:-1.2799   1st Qu.:-1.2499   1st Qu.:-1.2283  
##  Median : 0.5684   Median : 0.4237   Median : 0.5033   Median : 0.4981  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.7789   3rd Qu.: 0.7504   3rd Qu.: 0.7518   3rd Qu.: 0.7396  
##  Max.   : 1.7260   Max.   : 2.7328   Max.   : 2.3945   Max.   : 3.5285

Todas las medias son 0, veamos las desviaciones estandard por columna:

apply(data.norm, 2, sd)

##   Wing Weight Culmen Hallux 
##      1      1      1      1

data.norm.df = as.data.frame(data.norm)

ggpairs(data.norm.df, upper = list(continuous = "density")) +
  theme_bw()

De forma general podemos ver que ninguna variable de interés tiene una distribución normal, son distribuciones bimodales con media 0 y sd 1.

1.4 Modelado

Búsqueda de clústeres para k-means

La desventaja de los métodos de codo y silueta promedio es que solo miden una característica de agrupación global. Un método más sofisticado es usar la estadística de brecha que proporciona un procedimiento estadístico para formalizar la heurística de codo/silueta para estimar el número óptimo de grupos.

# Elbow method
fviz_nbclust(data.norm, kmeans, method="wss") + geom_vline(xintercept=4, linetype=2)+ labs(subtitle = "Elbow method")

# Silhouette method
fviz_nbclust(data.norm, kmeans, method = "silhouette")+ labs(subtitle = "Silhouette method")

# Gap statistic
set.seed(123)
fviz_nbclust(data.norm, kmeans, nstart = 25,  method = "gap_stat", nboot = 10)+ labs(subtitle = "Gap statistic method")

Podemos ver cómo diferentes métodos nos dan 4 ó 2 clusters para Kmeans. La desventaja de los métodos de codo y silueta promedio es que solo miden una característica de agrupación global. Un método más sofisticado es usar la estadística de brecha que proporciona un procedimiento estadístico para formalizar la heurística de codo/silueta para estimar el número óptimo de grupos.

Vamos hacer un promedio con diferentes métodos, el paquete NbClust evalúa el promedio número de clústeres más apropiado.

nb <- NbClust(data.norm, distance = "euclidean", min.nc = 2,
              max.nc = 10, method = "kmeans")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
##

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 12 proposed 2 as the best number of clusters 
## * 3 proposed 3 as the best number of clusters 
## * 6 proposed 4 as the best number of clusters 
## * 1 proposed 6 as the best number of clusters 
## * 2 proposed 8 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  2 
##  
##  
## *******************************************************************

El número óptimo de clusters nos dió 2, vamos a calcular y comprobar los resultados del algoritmo con 2 clusters en tándem para cada variable entre datos normalizados y datos antes de normalizar.

Aplicación de k-means con 2 clusters

# clusters para los datos sin normalizar 
data_kmeans <- kmeans(data[2:5], 2)
# clusters para los datos normalizados 
data.norm_kmeans2 <- kmeans(data.norm.df, 2)

# par() to create multiple plots at once.
# mfrow() to create a multi-paneled plotting window. 
par(mfrow=c(1,2))

#Sacamos las gráficas para los datos sin normalizar
for (i in 1:3){
  for (j in (i+1):4){
      plot1 <- ggplot(data = data.norm.df, aes(x=get(colnames(data.norm.df)[i]), 
                                               y=get(colnames(data.norm.df)[j]), 
                                               color = data$Species))+
        geom_point()+
        ggtitle("Puntos data real")+
        xlab(colnames(data.norm)[i])+
        ylab(colnames(data.norm)[j])+
        theme(legend.position="bottom")
    
      plot2 <- ggplot(data = data.norm.df, aes(x=get(colnames(data.norm.df)[i]), 
                                               y=get(colnames(data.norm.df)[j]), 
                                               color = factor(data.norm_kmeans2$cluster)))+
      geom_point()+
      scale_colour_manual(values = c("red", "blue", "green","black"))+
      ggtitle("Puntos cluster K-means")+
      xlab(colnames(data.norm.df)[i])+
      ylab(colnames(data.norm.df)[j])+
      theme(legend.position="bottom")
    
    grid.arrange(plot2, plot1, nrow = 1)
  }
}

El grupo formado por los puntos negros de la clasificación K-means coincide con los puntos RT y SS de la clasificación real. Podríamos enumerar los grupos como

grupo 1 : color azul ( RT )
grupo 2 : color rojo (CH y SS)

validar el método de agrupamiento.

Después de ajustar los datos en conglomerados utilizando diferentes métodos de conglomerados, deseamos medir la precisión del conglomerado. En la mayoría de los casos, usamos métricas intraclúster o interclúster como medidas. Cuanto mayor sea la distancia entre grupos, mejor será, y cuanto menor sea la distancia entre grupos, mejor será. La medida within.cluster.ss representa la suma de cuadrados dentro de los grupos y avg.silwidthrepresenta el ancho promedio de la silueta .

within.cluster.ss la medición muestra qué tan estrechamente relacionados están los objetos en grupos; cuanto menor sea el valor, más objetos estrechamente relacionados estarán dentro del clúster.
avg.silwidth una medida que considera qué tan estrechamente relacionados están los objetos dentro del grupo y cómo se separan los grupos entre sí. El valor de la silueta suele oscilar entre 0 y 1; un valor más cercano a 1 sugiere que los datos están mejor agrupados.

Veamos puntuaciones para data.norm_kmeans2:

# cluster.stats : Computes a number of distance based statistics, which can be used for cluster validation,
# avg.silwidth : silhouette information according to a given clustering in  clusters.

score = cluster.stats(dist(data.norm), data.norm_kmeans2$cluster)
score[c("within.cluster.ss","avg.silwidth")]

## $within.cluster.ss
## [1] 486.3121
## 
## $avg.silwidth
## [1] 0.7764635

Aplicación k-means 3 cluster

# clusters para los datos normalizados 
data.norm_kmeans3 <- kmeans(data.norm.df, 3)

# par() to create multiple plots at once.
# mfrow() to create a multi-paneled plotting window. 
par(mfrow=c(1,2))

#Sacamos las gráficas para los datos sin normalizar
for (i in 1:3){
  for (j in (i+1):4){
      plot1 <- ggplot(data = data.norm.df, aes(x=get(colnames(data.norm.df)[i]), 
                                               y=get(colnames(data.norm.df)[j]), 
                                               color = data$Species))+
        geom_point()+
        ggtitle("Puntos data real")+
        xlab(colnames(data.norm.df)[i])+
        ylab(colnames(data.norm.df)[j])+
        theme(legend.position="bottom")
    
      plot2 <- ggplot(data = data.norm.df, aes(x=get(colnames(data.norm.df)[i]), 
                                               y=get(colnames(data.norm.df)[j]), 
                                               color = factor(data.norm_kmeans3$cluster)))+
        geom_point()+
        scale_colour_manual(values = c("red", "blue", "green","black"))+
        ggtitle("Puntos cluster K-means")+
        xlab(colnames(data.norm.df)[i])+
        ylab(colnames(data.norm.df)[j])+
        theme(legend.position="bottom")
    
      grid.arrange(plot2, plot1, nrow = 1)
  }
}

Hemos creado los gráficos con los datos nomalizados para k-means, podemos ver como los puntos no están bien delimitados en referencia con los datos reales. Veamos puntuaciones para data.norm_kmeans3:

score = cluster.stats(dist(data.norm), data.norm_kmeans3$cluster)
score[c("within.cluster.ss","avg.silwidth")]

## $within.cluster.ss
## [1] 389.0211
## 
## $avg.silwidth
## [1] 0.5021049

2

Con el juego de datos proporcionado se realiza un estudio DBSCAN y OPTICS

minPts define la mínima densidad aceptada alrededor de un centroide.
e (eps) es el radio de nuestras vecindades alrededor de un punto de datos p .

# Distance matrix
dd <- daisy(data.norm)
str(dd)

##  'dissimilarity' num [1:405450] 0.132 0.255 2.031 3.635 0.657 ...
##  - attr(*, "Labels")= chr [1:901] "1" "2" "3" "4" ...
##  - attr(*, "Size")= int 901
##  - attr(*, "Metric")= chr "euclidean"

quantile(dd)

##         0%        25%        50%        75%       100% 
## 0.01052301 0.62847090 1.78228063 3.92520132 7.41105837

2.1 OPTICS

Calculamos OPtics para vecindad 10 y eps=0.8

res <- optics(data.norm, minPts = 10, eps = 1)
res

## OPTICS ordering/clustering for 901 objects.
## Parameters: minPts = 10, eps = 1, eps_cl = NA, xi = NA
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
##                   eps_cl, xi

Gráfica de alcanzabilidad

plot(res)

Sabemos que tenemos tres especies de halcones.

Determinación del valor óptimo de e (eps): El método propuesto aquí consiste en calcular las distancias de los k-vecinos más cercanos en una matriz de puntos. El valor de k corresponde a MinPts . A continuación, k -distancias se trazan en orden ascendente y el objetivo es determinar “la rodillael”codo” que corresponde al parámetro e óptimo.

Un codo corresponde a un umbral en el que se produce un cambio brusco a lo largo de la curva de distancia k . La función kNNdistplot() se puede utilizar para dibujar la gráfica de distancia k .

# Fuente : https://livebook.manning.com/book/machine-learning-with-r-the-tidyverse-and-mlr
kNNdistplot(data.norm, k = 10)
abline(h = c(0.4, 0.8))

Esta región donde la curva se inclina hacia arriba es el valor óptimo de épsilon (codo) a la distancia vecina más cercana en esta inflexión. Cuando la curva augmenta nos estamos alejando de la densidad y entremos es los valores de baja densidad que son los valores outlaiers. Usando este método, seleccionamos 0.4 y 0.8 como los límites inferior y superior sobre los cuales afinar el epsilon.

2.2 DBSCAN

Establezco el parámetro eps en eps_cl = 0.4.

res <- extractDBSCAN(res, eps_cl = .4)
res

## OPTICS ordering/clustering for 901 objects.
## Parameters: minPts = 10, eps = 1, eps_cl = 0.4, xi = NA
## The clustering contains 3 cluster(s) and 56 noise points.
## 
##   0   1   2   3 
##  56 544  56 245 
## 
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
##                   eps_cl, xi, cluster

Podemos ver como para eps_cl = 0.4 tenemos tres clusters y 56 puntos de ruido. Veamos el plot

plot(res)

Observamos en el gráfico anterior como se han coloreado los 3 clusters, en negro se mantienen los valores extremos.

hullplot(data.norm, res)

Podemos observar la cantidad de puntos que no se han clasificado.

Veamos ahora la validación del método:

score = cluster.stats(dist(data.norm), res$cluster)
score[c("within.cluster.ss","avg.silwidth")]

## $within.cluster.ss
## [1] 448.3559
## 
## $avg.silwidth
## [1] 0.6476835

Vamos a provar de clasificar con diferentes valores de eps. Probamos para eps_cl=0.8

res <- extractDBSCAN(res, eps_cl = .8)
res

## OPTICS ordering/clustering for 901 objects.
## Parameters: minPts = 10, eps = 1, eps_cl = 0.8, xi = NA
## The clustering contains 2 cluster(s) and 26 noise points.
## 
##   0   1   2 
##  26 560 315 
## 
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
##                   eps_cl, xi, cluster

Podemos observar como el algoritmo ha producido 2 clusters con más los 20 outliers.

plot(res)

score = cluster.stats(dist(data.norm), res$cluster)
score[c("within.cluster.ss","avg.silwidth")]

## $within.cluster.ss
## [1] 495.0134
## 
## $avg.silwidth
## [1] 0.7018825

Ahora, habiendo analizado diferentes medidas para eps_cl y sabiendo que el data set original cuenta con 3 especies de haclones, nos aseguramos los clusters con mínimos outliers, graficamos para eps_cl=0.6 ya que genera 3 clusters con los outliers más bajos.

res <- extractDBSCAN(res, eps_cl=0.6)
res

## OPTICS ordering/clustering for 901 objects.
## Parameters: minPts = 10, eps = 1, eps_cl = 0.6, xi = NA
## The clustering contains 3 cluster(s) and 36 noise points.
## 
##   0   1   2   3 
##  36 556  60 249 
## 
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
##                   eps_cl, xi, cluster

Tenemos puntos 36 aouliers, un cluster con 556 puntos, un segundo grupo con 60 y un tercer grupo con 249 puntos.

plot(res)

hullplot(data.norm, res)

Como hemos hecho anteriormente, pasamos a comprobar como se comporta el algoritmo en tándem para cada variable entre datos normalizados y datos reales.

for (i in 1:3){
  for (j in (i+1):4){
    plot1 <- ggplot(data = data.norm.df, aes(x=get(colnames(data.norm)[i]), 
                                       y=get(colnames(data.norm)[j]), color = data$Species))+
      geom_point()+
      ggtitle("Valores reales")+
      xlab(colnames(data.norm)[i])+
      ylab(colnames(data.norm)[j])
    
    plot2 <- ggplot(data = data.norm.df, aes(x=get(colnames(data.norm.df)[i]), 
                                       y=get(colnames(data.norm.df)[j]), color = factor(res$cluster)))+
      geom_point()+
      scale_colour_manual(values = c("red", "blue", "green","black"))+
      ggtitle("Valores del cluster")+
      xlab(colnames(data.norm.df)[i])+
      ylab(colnames(data.norm.df)[j])
    
    grid.arrange(plot2, plot1, nrow = 1)
  }
}

Con el procedimiento DBSCAN los puntos se reproducen mejor que con K-means.

Grupo 0 : valores atípicos
Grupo 1 : valores de especie TR
Grupo 2 : valores de especie CH
Grupo 3 : valores de especie SS

score = cluster.stats(dist(data.norm), res$cluster)
score[c("within.cluster.ss","avg.silwidth")]

## $within.cluster.ss
## [1] 380.9451
## 
## $avg.silwidth
## [1] 0.6740071

3

3.1 Comparación resultados K-means y DBSCAN

Probamos el algoritmo k-means para 2 y 3 clusters; medimos el resultado en ambos casos y para 2 clusters la medición de calidad de agrupamiento fue de 0.776 y para 3 clusters fue de 0.502 o sea una medida de agrupación bastante por debajo para 3 clusters que ya sabemos que originalmente así sería.

K-means 3 clusters: 0.502
K-means 2 clusters: 0.776

Posteriormente hicimos varios análisis con DBSCAN.

DBSCAN eps_cl=0.4 : 0.647 con 3 cluster
DBSCAN eps_cl=0.8 : 0.708 con 2 cluster
DBSCAN eps_cl=0.6 : 0.674 con 3 cluster

Podemos ver que el que tuvo mejor puntación de calidad de agrupamiento fue para kemans con dos clusters. Si comparamos para tres clusters K-means queda con peor score siendo DBSCAN eps_cl=0.6 : 0.674 con 3 cluster la mejor opción.

3.2 Pros y contras de K-means y DBSCAN

En el algoritmo k-means se necesita predeterminar el número de clusters; esto no sucede con DBSCAN ya que se obtiene mediante la baja o alta densidad de puntos concentrados en áreas.

Para DBSCAN necesitamos determinar la densidad apropiada para los parámetros epsilon y minpts, esta tarea es más compleja porque depende del conocimiento del dominio y el conjunto de datos que se aportan al algoritmo.

El algortimo k-means crea clusters de forma más o menos esférica o convexa y deberían tener el mimso tamaño de caracteristicas, en cambio DBSCAN los clusyters que forman son más aleatorios y sus tamaños no tiene por qué ser igual a sus caracteristicas, odea es capaz de identificar clusteres de cualquier forma geométrica, en cambio k-means sólo lo hace de forma circular.

3.3 Conclusiones

Comenzamos el trabajo con un conjunto de datos proporcionado en la PEC con 908 obs. of 19 variables. donde se miden de forma aleatoria diferentes atributis de tres especies de halcones. Hemos trabajado métodos de clasificación no supervisada para cuatro variables numéricas (Wing, Weight, Culmen y Hallux) sabiendo de ante mano cuál era la clasificación por especie. Por cada método de clasificación trabajado se ha evaluado su calidad de agrupamiento tanto de centro como de lejanía. Estas medidas nos han proporcionado información para saber qué tan bien o mal se clasifican y basándonos en esto la mejor clasificación obtenida ha sido para k-mean con dos clusters con un puntaje avg siluette de 0.7764635, la decisión para evaluar en dos clusters deriva del análisis preliminar de consenso. La sengunda mejor puntuación avg siluette es para DBSCAN eps_cl=0.6 : 0.674 con 3 cluster. K-means para 3 clusters ha proporcinado un score bastante más bajo de 0.502

4

4.1 Ventajas k-means y DBSCAN

K-MEANS :

VENTAJAS

Relativamente fácil de implementar.
Escala a grandes conjuntos de datos.
Garantiza la convergencia.
Puede iniciar en caliente las posiciones de los centroides.
Se adapta con facilidad a los ejemplos nuevos.
Se generaliza a clústeres de diferentes formas y tamaños, como clústeres elípticos.

DESVENTAJAS

Elige K clusters manualmente.
Depende de los valores iniciales.
No es bueno con el agrupamiento de datos de diferentes tamaños y densidades.
Valores atípicos del agrupamiento en clústeres.
Escalamiento con cantidad de dimensiones

APLICACIONES

Geostático
Visión computacional
Segmentación de mercado Estudio de terremotos
Uso de suelo

DBSCAN :